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A  complete,  unified  description  is  given  of  the  design,  implementation  and  use  of  a  family 
of  very  fast  and  efficient  large-scale  minimum-cost  (primal  simplex)  network  programs.  The 
class  of  capacitated  generalized  transshipment  problems  solved  includes  the  capacitated  and 
uncapacitated  generalized  transportation  problems  and  the  continuous  generalized  assignment 
problem,  as  well  as  the  pure  network  flow  models  which  are  specializations  of  these  problems. 
These  formulations  are  used  for  a  large  number  of  diverse  applications  to  determine  how  (or 
at  what  rate)  flows  through  the  arcs  of  a  network  can  minimize  total  shipment  costs.  A 
generalized  network  problem  can  also  be  viewed  as  a  linear  program  with  at  most  two  nonzero 
entries  in  each  column  of  the  constraint  matrix;  this  property  is  exploited  in  the  mathematical 
presentation  with  special  emphasis  on  data  structures  for  basis  representation,  basis  manipula¬ 
tion,  and  pricing  mechanisms.  A  literature  review  accompanies  computational  testing  of 
promising  ideas,  and  extensive  experimentation  is  reported  which  has  produced  GENNET,  an 
extremely  efficient  family  of  generalized  network  systems. 

(NETWORKS/GRAPHS— FLOW  ALGORITHMS;  NETWORKS/GRAPHS— GEN¬ 
ERALIZED) 


1.  Introduction 


This  paper  reports  the  development  of  a  large-scale  primal  network  code  for  solving  capacitated 
generalized  transshipment  problems.  The  capacitated  generalized  transshipment  problem  is  the  most  general 
of  the  minimum  cost  flow  models  in  continuous  variables,  which  include  the  capacitated  and  uncapacitated 
transportation  problems  and  the  continuous  generalized  assignment  problem  as  well  as  the  pure  network 
specializations  of  these  problems.  These  models  are  used  for  a  large  number  of  diverse  applications  that 
include  transportation  of  goods,  design  of  reservoir,  communications,  and  pipeline  systems,  assignment  of 
personnel  and  machinery  to  jobs,  bid  evaluation,  currency  exchange  and  cash  management,  production, 
sales  and  inventory  planning,  and  many  others.  For  further  discussion  of  these  applications  see  survey 
articles  such  as  Bradley  (1975),  or  Glover  et  al.  (1977,  1978),  and  textbooks  (as  well  as  their  cited  references) 
such  as  Dantzig  (1963),  Jensen  and  Barnes  (1980),  and  Kennington  and  Helgason  (1980). 

The  capacitated  generalized  transshipment  model  and  its  specializations  are  minimum-cost  network  flow 
problems.  The  goal  is  to  determine  how  (or  at  what  rate)  flows  through  the  arcs  of  a  network  can  minimize 
shipment  costs.  The  network  is  a  directed  graph  G,  defined  by  a  set  of  nodes,  N,  and  a  set  of  arcs  A,  with 
ordered  pairs  of  nodes  (tail,  head)  as  elements  indexed  by  k  :  (fk,  tk).  For  each  arc  there  is  a  shipping  cost 
per  unit  flow,  ck ,  a  minimum  allowable  flow  (or  lower  bound),  lk ,  and  a  maximum  allowable  flow  (or  upper 
bound,  or  capacity),  uk.  In  addition,  there  are  coefficients  (or  multipliers,  or  gains,  or  losses),  ak  and  bk 
which  can  change  the  magnitude  of  (or  amplify,  or  attenuate)  each  unit  of  flow  respectively  entering  and 
leaving  arc  k.  Each  node  is  either  a  supply  node  where  units  of  the  good  enter  the  network,  a  demand  node 
where  units  leave,  or  a  transshipment  node.  The  problem  is  to  minimize  total  costs  with  flows,  xk,  that 
satisfy  the  associated  lower  and  upper  bounds  and  preserve  the  conservation  of  flow  at  each  node: 

MIN  2  ck*k  (GN) 

kE  A 


where: 


s.t.  2  £***+  S  ieN, 

kE  A_  kE  A_ 

with  fk  =  i  with  tk  =  i 

4  <  xk  <  uk  ,  k  e  A_ , 


b,  = 


supply  if  i  is  a  supply  node; 

-  demand  if  i  is  a  demand  node; 
0  otherwise. 


•Accepted  by  Donald  Erlenkotter,  former  Departmental  Editor;  received  April  21,  1982.  This  paper  has 
been  with  the  authors  5  {  months  for  2  revisions. 
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Note  that  any  linear  program  with  at  most  two  nonzero  coefficients  associated  with  each  variable  is  a 
generalized  network  (GN).  Formulation  (GN)  is  a  generalized  transshipment  model  (GT)  if  all  are  +  1,  in 
which  case  the  corresponding  coefficient  is  called  the  multiplier  mk  =  -  bk .  For  purposes  of  exposition,  we 
will  address  (GT)  since,  assuming  the  existence  of  a  finite  upper  bound  on  each  variable,  it  is  possible  to 
transform  the  (GN)  coefficients — by  scaling  or  by  reflecting  variables  with  respect  to  their  upper  bounds — so 
that  one  coefficient  is  + 1  for  each  variable.  The  other  coefficient  for  each  variable  then  becomes  the 
multiplier  value,  mk,  and  the  generalized  transshipment  (GT)  network  flow  interpretation  results  with  a  node 
for  each  constraint  and  a  directed  arc  for  each  variable.  If  a  variable  has  two  nonzero  coefficients,  its  arc  is 
directed  away  from  the  node  corresponding  to  the  constraint  with  the  +  1  coefficient;  a  variable  with  just 
one  nonzero  coefficient  ( ±  1 )  corresponds  to  an  arc  forming  a  self-cycle,  leaving  and  returning  to  the  same 
node. 

Some  generalized  networks  (GN),  such  as  those  obtained  by  relaxing  integer  restrictions  on  flow  variables, 
cannot  be  scaled  conveniently  to  (GT).  These  (GN)  are  accommodated  by  obvious  minor  modifications  in 
the  following  (GT)  presentation.  We  have  developed  codes  for  solving  (GN)  and  specialized  them  for  (GT) 
problems. 

Generalized  networks  can  be  solved  as  linear  programming  problems,  but  contemporary  commercial 
linear  programming  systems  consume  much  more  computer  time  and  data  storage  region  than  special 
purpose  network  codes.  Indeed,  the  advantage  of  network  codes  is  so  pronounced  that  it  is  even  worthwhile 
to  develop  special  linear  programming  procedures  to  exploit  intrinsic  network  structure  found  embedded 
within  more  general  models  (e.g.,  Kennington  and  Helgason  1980,  McBride  1981,  Brown  and  Graves  1975, 
and  Brown  and  McBride  1981).  Brown  and  Wright  (to  appear)  and  Brown,  McBride  and  Wood  (1983)  show 
that  many  real-life  linear  programs  contain  a  large  embedded  generalized  network  structure. 

These  models  are  widely  used  because  they  accurately  describe  a  large  variety  of  important  applications. 
Generalized  networks  not  only  directly  represent  gains  with  mk  >  1  (e.g.  interest  return  on  investment,  heat 
gain,  etc.)  and  losses  with  0  <mk  <  1  (e.g.,  evaporation,  voltage  drop,  attrition,  etc.)  in  the  flows,  but  also 
admit  conversions  of  units  for  these  flows  (e.g.,  machine  time  to  output  pieces,  lira  to  yen,  etc.).  In  addition, 
mk  <  0  represents  situations  without  obvious  physical  flow  interpretation  (e.g.,  flows  which  “enter  the  head” 
of  arc  k  in  proportion  mk  of  those  entering  the  tail),  but  which  nonetheless  provide  valuable  modelling  tools. 
There  has  been  continuing  growth  of  interest  in  network  models  because  efficient  computer  programs  have 
made  possible  the  reliable,  economic  solution  of  problems  with  more  variables  than  virtually  any  other 
optimization  technique  (e.g.,  the  pure  network  system,  GNET  (Bradley  et  al.  1977),  has  been  installed  at 
hundreds  of  sites  worldwide  and  is  now  cited  as  a  routine  research  tool).  Perhaps  most  important,  networks 
are  readily  accepted  by  nonanalysts  and  are  consequently  extremely  popular  operations  research  models. 

Although  several  papers  have  been  written  in  this  general  area,  and  significant  computational 
breakthroughs  have  been  reported,  there  has  not  previously  been  a  single,  unified  description  of  a  complete 
implementation,  nor  have  “new  generation”  computer  programs  been  made  generally  available  to  the 
academic  community.  Here  we  report  the  research  and  computational  experiments  which  have  produced 
GENNET,  an  extremely  efficient  family  of  network  optimization  systems.  GENNET  exploits  pure  network 
structure  embedded  in  generalized  networks,  and  specializes  intrinsically  to  GNET  (Bradley  et  al.  1977), 
when  both  systems  are  applied  to  pure  networks,  (the  floating  point  arithmetic  in)  GENNET  requires  about 
15  percent  more  time  than  GNET.  An  important  objective  of  this  paper  is  to  make  these  new  approaches 
easily  accessible  to  a  wide  audience  via  a  clear  exposition  and  concrete  examples  of  efficient  FORTRAN 
programs.  Further,  the  availability  of  the  (GT)  and  (GN)  computer  programs  will  now  make  it  possible  for 
other  investigators  to  reproduce  and  extend  our  experimental  results. 

Bradley,  Brown  and  Graves  (1977)  trace  the  historical  developments  leading  to  contemporary  primal 
simplex  pure  network  algorithms,  their  supporting  data  structures,  and  efficient  implementations  such  as 
GNET.  For  other  subclasses  of  generalized  networks,  algorithms  have  been  reported  by  Jewell  (1962), 
Eisemann  (1964),  Maurras  (1972),  Glover,  Hultz,  Klingman  and  Stutz  (1977,  1978,  1974),  Balachandran 
(1976),  and  Jensen  and  Bhaumik  (1977).  An  efficient  algorithm  for  large  generalized  network  problems  has 
been  developed  by  Glover,  Klingman,  Hultz,  Stutz,  Karney  and  Elam  (1979,  1977,  1978,  1973,  1973). 
However,  their  contributions  are  scattered  among  the  papers  referenced;  Kennington  and  Helgason  (1980) 
and  Jensen  and  Barnes  (1980)  provide  textbook  descriptions  of  computations  apparently  gleaned  from  these 
papers,  providing  an  adequate  treatment  of  the  graphical  algorithm  with  more  computational  advice  than 
the  seminal  presentation  by  Dantzig  (1963)  but  few  details  of  efficient  basis  updating. 


2.  The  Approach 

Our  approach  continues  with  generalized  networks  in  precisely  the  philosophical 
vein  of  the  pure  network  exposition  of  Bradley,  Brown  and  Graves  (1977,  p.  3ff.):  we 
seek  data  structures  and  algorithms  that  yield  efficient  implementations  without 
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abandoning  the  flexibility  of  a  general  large-scale  mathematical  programming  perspec¬ 
tive.  We  introduce  few  of  the  details  of  the  general  bounded-variable  simplex  algo¬ 
rithm,  and  we  repeat  little  of  the  underlying  pure  network  material;  the  assiduous 
reader  might  well  review  the  prior  paper  for  which  this  is  an  intimate  companion. 

We  continue  with  a  brief  description  of  the  algebraic  specialization  of  the  simplex 
method  for  generalized  networks.  Specific  design  decisions  and  experiments  carried 
out  with  GENNET  are  described,  including  computational  tests  of  alternate  ap¬ 
proaches.  Some  extensions  of  GENNET  are  presented  to  further  exploit  special 
problem  structure. 

3.  Generalized  Network  Specialization 

Efficient  primal  simplex  specialization  to  the  generalized  network  case  depends  upon 
the  well-known  result  that  any  generalized  network  basis  can  be  put  in  nearly  (upper) 
triangular  form  by  simple  permutation  of  rows  and  columns.  This  inherent  near 
triangularity  can  be  exploited  by  direct  solution  of  the  simplex  equations  with  modified 
forward,  or  back  substitution.  Fortuitously,  this  basis  structure  also  leads  to  extremely 
fast  solution  updates  orchestrated  in  concert  with  efficient  dynamic  reorganization  of 
each  new  basis. 

Theorem  (e.g.,  Dantzig  1963).  Any  basis  B  extracted  from  a  generalized  network 
problem  can  be  put  in  the  form  (1 )  by  rearranging  rows  and  columns: 

B 1 

B1 

B‘ 

BP 

where  each  square  submatrix  component  Bl  is  either  upper  triangular  or  nearly  upper 
triangular  with  only  one  element  below  the  diagonal. 

Consider  the  generalized  network  problem  given  in  Figure  1,  a  generalized  transship¬ 
ment  problem  with  2  sources,  10  sinks,  15  nodes,  and  30  arcs.  We  will  use  this  problem 
to  illustrate  concepts  and  efficient  solution  methods.  In  this  problem  the  multipliers  are 
all  positive.  For  arc  k  directed  from  node  i  to  node  j,  if  flow  xk  leaves  node  i,  then 
mkxk  arrives  at  node  j.  When  mk  >  1  the  amount  arriving  at  node  j  is  greater  than  the 
amount  leaving  node  i.  This  would  be  the  case  in  cash  flow  problems  when  the  arc 
corresponds  to  an  investment  and  mk  =  (1  +rk)  with  rk  the  rate  of  return.  When 
mk  <  1  then  the  amount  arriving  at  node  j  is  less  than  the  amount  leaving  node  /.  Here 
the  loss  could  correspond  to  evaporation,  taxation,  transmission  loss,  brokerage  fees, 
seepage  or  deterioration.  Pure  network  arcs  are  indicated  by  mk  =  1  (and  may  be  even 
more  abundant  in  real  problems  than  in  our  example).  For  clarity,  minimum  allowable 
flows  are  zero,  and  maximum  allowable  flows  are  not  specified. 

Figure  2  shows  a  basis  for  the  problem  introduced  in  Figure  1.  (In  this  simple  case, 
there  is  only  one  component:  p  =  1.) 

A  unique  subgraph  partition  of  G  denoted  GB  corresponds  to  B.  Let  AH  =  { a  7 1  a  J  is 
an  arc  associated  with  bj,  a  column  of  B),  then  GB  =  [K,Ab\  denotes  the  directed 
graph  associated  with  the  basis  B.  To  each  submatrix  Bl  of  B  there  corresponds  a 
component  of  GB  denoted  by  G1  =  [A^, A1].  N1  is  the  set  of  nodes  corresponding  to  the 
rows  of  Bl  and  A1  is  the  set  of  arcs  corresponding  to  the  columns  of  Bl.  It  is  known 
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Figure  1.  A  Single  Commodity  Generalized  Transshipment  Problem  (GT). 


(e.g.,  Dantzig  1963)  that  G1  is  either  a  rooted  tree  or  a  one-tree  (a  tree  with  an 
additional  arc  forming  one  cycle).  If  Gl  is  a  rooted  tree,  then  B 1  is  upper  triangular.  If 
G1  is  a  one-tree,  then  Bl  is  upper  triangular  except  for  one  element  below  the  diagonal. 
Each  component  can  be  viewed  as  having  one  cycle  if  we  assert  that  each  rooted  tree 
has  a  self-cycle  corresponding  to  its  root  node. 

In  our  example  basis,  the  subgraph  GB  has  only  one  component,  (one-tree)  shown  in 
Figure  2:  GB  is  a  one- tree,  with  nodes  2,  3  and  11  composing  its  cycle. 

As  in  the  pure  network  case,  this  near  triangulation  and  associated  subgraph  are 
naturally  represented  by  a  predecessor  function  p(  ),  and  a  predecessor  graph  (which 
does  not  preserve  the  orientation  of  arcs  in  the  original  network).  The  predecessor 
function  can  be  used  iteratively  to  construct  the  unique  backpath  from  any  node  to  the 
root  (or  cycle);  the  backpath  includes  all  nodes  on  the  cycle.  The  immediate  successors 
of  a  node,  if  any,  are  the  first  nodes  encountered  on  all  paths  except  the  backpath  to 
the  root,  and  all  the  nodes  on  these  paths  are  called  successors. 

Note  that  each  basis  may  have  many  near  triangulations.  However,  all  such  near 
triangulations  yield  the  same  predecessor  function  and  graph  (where  the  right  to  left 
ordering  of  successors  of  any  node  is  immaterial).  Thus,  the  predecessor  graph  does 
not  completely  represent  a  near  triangulation  without  additional  information:  an 
ordering  of  the  rows  (nodes).  For  algebraic  reasons,  we  restrict  such  partial  ordering  to 
preorder  (Bradley  et  al.  1977),  in  which  a  node  i  always  precedes  its  successors,  if  any, 
and  in  which  all  its  successors,  if  any,  precede  any  node  which  does  not  precede  node  i. 
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4.  Implementation 

For  didactic  reasons,  we  begin  by  introducing  a  complete  primal  generalized 
network  algorithm  using  a  preorder  traversal  method.  Controversial  alternatives  are 
deferred  until  this  paradigm  is  presented.  Hereafter,  notation  with  upper  case  italic 
letters  followed  by  parentheses  indicates  a  program  data  array.  For  instance,  the 
predecessor  array  is  referred  to  as  P( ). 

Static  arc  storage  is  used  for  tails  T( ),  head  H( ),  costs  C( ),  multipliers  MUL( ),  and 
capacities  CP( ).  Contiguous  storage  by  tail,  or  by  head  node  reduces  T( ),  or  H( )  to  a 
hierarchical  node-length  entry  point  array.  (GENNET  uses  contiguous  storage  by  head 
node,  as  does  GNET.)  Lower  bounds  on  arc  flow  are  translated  out  prior  to  solution, 
with  appropriate  adjustment  of  the  initial  right-hand  side  of  (GT),  and  of  CP(  ).  The 
sign  bit  of  CP( )  is  available  to  indicate  arcs  nonbasic  at  their  upper  bounds  ( reflected 
with  flow  —  CP( )). 

The  predecessor  function  and  its  array  P(  )  are  defined  so  that  the  basic  arcs  in  a 
cycle  are  oriented  uniformly  in  a  directed  cycle.  All  basic  arcs  not  on  a  cycle  are 
oriented  so  that  a  backpath  is  created  to  a  cycle.  To  obtain  this  orientation  the 
direction  of  some  arcs  must  be  reversed,  and  the  sign  bit  of  the  predecessor  array  is 
used  to  indicate:  if  P(I)  <  0,  then  the  orientation  of  arc  (I,  —  P(I))  is  reversed  from  its 
original  orientation.1 

A  depth  array  D(  )  reveals  for  each  node  the  number  of  nodes  on  the  backpath 
before  encountering  a  cycle.  Nodes  on  a  cycle  have  depth  zero.  Number  of  successors, 
or  preorder  distance  are  acceptable  substitutes  for  depth  (Bradley  et  al.  1977),  but  are 
not  discussed  here. 

A  preorder  traversal  array  IT(  )  is  maintained  so  that  all  preorder  successors  of  a 
cycle  node  are  encountered  before  another  cycle  node.  It  is  convenient  to  make  this  a 
circular  list  for  each  near  triangular  basis  component  by  setting  IT(  )  of  the  last 
preorder  node  in  the  component  equal  to  the  first  preorder  node  in  the  component. 

The  components  of  GB  are  not  interconnected,  or  equivalently,  the  submatrices  Bl 
in  (1)  do  not  have  common  rows  or  columns.  Consequently,  the  p  components  of  a 
basis  may  be  represented  in  a  single  set  of  node-length  arrays. 

The  array  X(  )  contains  the  values  of  basic  variables,  values  of  dual  variables  (or 
simplex  multipliers,  or  node  potentials)  are  stored  in  U(  ),  and  IVAR(  )  gives  the 


Node 

Predecessor 

Depth 

Traversal 

Basic  Variable 

Dual 

Cycle  Factor 

P() 

D() 

IT() 

IVAR(  ) 

X() 

U() 

FAC() 

1 

13 

2 

6 

24 

22.860 

16.665 

* 

2 

3 

0 

5 

2 

118.169 

47.148 

-0.48101 

3 

11 

0 

10 

21 

45.826 

31.678 

-0.48101 

4 

9 

2 

2 

13 

0.0 

5.418 

* 

5 

-2 

1 

8 

4 

0.0 

27.552 

* 

6 

-3 

1 

15 

7 

21.308 

-  3.782 

* 

7 

-2 

1 

11 

10 

3.434 

-  38.898 

* 

8 

-  2 

1 

14 

12 

27.087 

27.487 

* 

9 

—  3 

1 

4 

14 

9.380 

1.958 

* 

10 

-3 

1 

13 

17 

13.150 

28.606 

* 

11 

-2 

0 

3 

20 

4.170 

-  16.053 

-0.48101 

12 

-2 

I 

7 

23 

2.076 

32.431 

* 

13 

-3 

1 

1 

26 

11.956 

-  13.922 

* 

14 

-2 

1 

12 

29 

22.204 

29.580 

* 

15 

-  3 

1 

9 

30 

16.550 

-  35.942 

* 

Figure  3.  GENNET  Basis  Representation  Arrays  (for  Basis  in  Figure  2). 


'This  is  the  complement  of  the  discipline  used  in  GNET  (Bradley  et  al.  1977). 
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location  of  basic  variables  in  the  arc  arrays.  The  array  FAC(  )  contains  the  cycle 
factors,  defined  later.  Figure  3  shows  these  arrays  for  the  basis  given  in  Figure  2. 

Generalized  networks  do  not  exhibit  totally  unimodular  bases.  Consequently,  float¬ 
ing  point  representation  is  required  for  X(  ),  U(  )  and  FAQ  ),  and  is  desirable  for 
arc-length  arrays  Q  ),  MUL(  ),  and  CP(  ). 

Step  SI,  Priceout 

The  reduced  cost  for  nonbasic  arc  k,  oriented  from  fk  to  tk  is  (given  the  current  dual 
solution  u  and  column  Nk ): 

rk  =  ck-  uNk  =  ck  -  ufk  +  mku,k 
=  C(k)  —  U(fk )  +  MUL(k)*U(4). 

(If  CP(k)  <  0,  arc  k  is  reflected  and  the  sign  of  rk  is  reversed.)  At  most,  one 
multiplication,  addition  and  subtraction  are  required.  Note  that  the  multiplication  is 
unnecessary  if  \mk\  =  1;  further  specialization  is  possible  for  sets  of  priced  arcs  with 
common  attributes.  If  arc  A:  is  a  logical  arc  (slack,  artificial,  or  surplus  variable)  then 
C(k)  can  be  logically  generated,  rather  than  explicitly  stored,  and  rk  =  C(k)  ±  U(  fk), 
depending  upon  the  sign  of  P (fk). 

From  the  example, 

r21  =  C(27)  -  U(5)  +  MUL(27)*U(13) 

=  20.67  -  (27.552)  +  0.88(-  13.922) 

=  -19.133. 

This  variable  will  be  used  as  the  entering  nonbasic  variable  for  further  illustration. 
Step  S2,  Ratio  Test 

For  the  determination  of  the  arc  to  leave  the  basis,  the  system  of  equations 
BZk  =  Nk  must  be  solved  for  the  transformed  column  Zk.  (  —  Nk  is  used  if  arc  k  is 
reflected.)  Due  to  the  near  triangularity  of  B,  this  incoming  column  transformation  can 
be  combined  with  the  ratio  test  in  a  single  integrated  process. 

Suppose  that  Nk  has  two  nonzero  coefficients  representing  an  arc  oriented  from  fk  to 
4,  and  that  the  arc  is  not  reflected.  An  apparent  complication  arises  if  fk  and  4  are  in 
separate  components  of  B  in  (1),  say  Bs  and  B\  respectively.  In  this  case  two  disjoint 
subsystems  must  be  solved  and  the  results  added  to  determine  the  nonzero  elements  in 
Zk.  The  subsystems  are: 

BsQh  =  ef,  and  B‘Qh  =  -  mke,k 

(e,t  is  the  4th  unit  vector;  Qf‘  and  Q‘k  are  disjoint  components  of  Zk). 

This  complication  is  inconsequential.  In  order  to  see  this,  consider  solving  one  of  the 
systems:  B‘Q'k  =  —  mk <4 .  The  only  nonzero  elements  of  Qh  will  be  those  that 
correspond  to  the  nodes  in  the  backpath  from  node  4.  As  we  shall  see,  this  follows 
from  the  manner  in  which  the  coefficient  -mk  in  row  4  propagates  during  substitu¬ 
tion  solution. 

Suppose  that  G'  is  a  rooted  tree.  The  backpath  from  node  4  can  be  denoted  by 
iteration  of  the  predecessor  function:  tk,  p(tk),  p(p(tk)),  .  .  .  ,  root;  this  sequence  is 
shown  in  the  following  figure  as  d,d  —  1,  .  .  .  ,  0,  analogous  to  the  backpath  length 
remaining  to  the  root.  The  values  of  a  and  b  for  each  basic  arc  depend  upon  the 
original  orientation  of  the  arc,  given  by  the  sign  of  P(  ).  For  original  orientation, 
P(  )  >  0  implies  that  a  =  + 1  and  b  is  the  multiplier  value,  P(  )  <  0  reverses  these 
definitions. 
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The  triangular  system  corresponding  to  this  backpath  is: 


a0  *1 

0 

al 

0 

Q'k  = 

ad- 2  bd- i 

0 

ad- 1  bd 

0 

ad 

- vik 

0 

By  back  substitution,  its  solution  is: 


qd  = 


-mjc 

ad 


qs  = 


^S+rfs+l 

as 


for  S  =  d—  1,  .  .  .  ,  0. 


Now  suppose  that  Gl  is  a  one-tree.  Let  node  tk  be  on  the  cycle.  The  backpath  is 
tk,  p(tk),  p(p(tk)),  ■  ■  ■  ,  c,  with  p(c)  =  tk  and  length  p;  this  sequence  is  shown  in  the 
figure  as  —  1,  —2,  ...  ,  —  p,  analogous  to  the  backpath  length  beyond  the  cycle  start. 


P<P( j) ) 

o  o« — — o 

b  .  o  _  b  _  o_b,  a, 
-/»♦!  -3  -2  -2  -1  -1 


The  corresponding  near  triangular  system  is: 


a — p  b _p+ 1 

0 

a-P+  i 

0 

Q"  = 

7 

7 

<3 

0 

a_2  6_, 

0 

b-p  «-i 

1 

IS 

0 
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By  modified  back  substitution,  its  solution  is: 


?-i  “ 


'  bg+  i<7sh 


for  8  =  —  2,  .  .  .  ,  —  p,  where 


-p 


/=i-  n 

8=  -  1 


The  cycle  factor  (or  loop  factor  Dantzig  1963),/,  is  common  to  all  nodes  in  the  cycle 
and  can  be  computed  when  the  cycle  is  created  and  stored  in  FAC( )  for  all  nodes  on 
the  cycle  so  that  it  is  immediately  available  at  this  step. 

Suppose  that  the  entering  arc  is  ( fk ,  tk)  with  coefficient  entries  ( ak ,  bk)  and  that  the  fk 
and  4  backpaths  converge.  Using  nomenclature  for  the  backpath  sequences  intro¬ 
duced  above,  the  corresponding  system  is: 


By  back  substitution,  its  solution  is: 


(  d  begins  tk  backpath  sequence), 
for  8  =  d,  d—  1,  .  .  .  ,  r, 

(d  begins  fk  backpath  sequence), 
for  8  =  d,  d  —  1,  .  .  ,  *  w, 


<li-\ 

% 


Suppose  that  entering  arc  (fk,tk)  with  coefficients  (ak,bk)  enters  and  that  the 
backpaths  converge  on  a  cycle.  The  corresponding  system  is: 


l+d- 
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By  modified  back  substitution,  its  solution  is: 

qd  =  bk/ a d  (d  begins  tk  backpath  sequence), 

~  bg+ 1?«+ 1 


qs  = 


for  8  =  d,  d—  1,  .  .  .  ,  r, 


-  b,qr  i 

—  1 

«-.v- I  / 

qs  = 

—  bs+ 1  qs+ 1 

<h 

qd  = 

ad 

qs  = 

~  bs+iqs+t 

a* 

for  8  =  —  s  —  2,  —  1 s, 

( d  begins  fk  backpath  sequence), 
for  8  =  d,  .  .  .  ,w, 


A 

i 

■is 

S 

i 

qi- 1  = 

A 

/’ 

qs  = 

bs+ 1  qs+ 1 

for 

as 

qs  = 

qs +  qs 

for 

For  the  cycle  in  Figure  2, 


8  =  -  l  -  2,  ...  ,  -  p,  -  1 ,  .  .  .  ,  -  s,  ...  ,  -  /, 

8=  -1,  ...,  —  p. 


/=  1- 


(-1) 


X 


(1) 


0.79 


X 


(—117) 

1 


0.48101. 


By  now  it  should  be  apparent  that  one  composite  back  substitution  scheme  will 
suffice  for  all  cases.  The  cycle  factor  is  applied  once  when,  and  if,  a  cycle  is 
encountered  on  a  backpath. 

If  the  backpaths  of  fk  and  tk  converge,  let  join  be  the  first  node  on  the  tk  backpath 
that  is  also  on  the  fk  backpath.  If  the  backpaths  converge  on  a  cycle  and  if  the  leaving 
arc  precedes  the  join  on  the  fk  backpath,  define  join  as  the  first  cycle  node  encountered 
on  the  fk  backpath.  If  the  backpaths  do  not  converge,  join  =  0. 

Several  schemes  are  available  for  identifying  the  join  efficiently  (Bradley  et  al.  1977). 
The  depth  (or  number  of  successors,  or  preorder  distance)  of  nodes  on  the  backpaths 
can  be  used  to  avoid  iterating  either  backpath  past  the  join.  Depth,  the  number  of 
nodes  on  the  backpath  until  a  cycle  is  encountered,  can  be  used  to  indicate  which 
backpath  node  is  deeper  and  should  be  iterated.  When  both  backpath  nodes  have 
matching  depths,  zero  depths  indicate  that  each  backpath  is  on  a  root  cycle.  By 
remembering  the  first  root  cycle  node  on  each  backpath,  further  iteration  will  either 
reveal  the  root  cycles  to  be  distinct  with  join  =  0,  or  coincident  with  join  defined  as 
above.  When  both  backpath  nodes  have  matching  depths  greater  than  zero,  the  nodes 
are  compared  for  equality.  A  match  indicates  the  join,  and  a  mismatch  indicates  that 
both  backpaths  should  be  iterated  for  another  comparison. 

If  a  join  is  encountered,  the  backpaths  have  merged  and  either  one  can  be  used  to 
complete  the  ratio  test  (GENNET  continues  the  tk  backpath).  When  a  cycle  is 
encountered  the  q  values  are  computed  on  the  first  pass  around  the  cycle.  On  the 
second  pass  around  the  cycle  the  q  values  are  computed  and  the  ratio  test  is 
completed. 

As  the  backpaths  are  iterated,  the  column  transformation  is  applied  and  the 
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resulting  terms  of  transformed  column,  Zk,  are  used  in  the  ratio  test,  seeking  the 
minimum  ratio : 


CP(/fc) 

X(0 


MIN- 


CP(IVAR(/))-X(/) 


the  capacity  of  the  incoming  arc, 
for  z,k  >  0,  node  /, 

for  Z/k  <  0. 


If  a  zero  ratio  is  encountered  during  this  process,  the  ratio  test  may  be  preemptively 
terminated. 


Step  S3,  Pivot 

If  CP(&:)  is  selected  as  the  minimum  ratio,  then  the  entering  variable  remains 
nonbasic  and  is  reflected  to  its  opposite  bound.  Only  the  flows  X( )  need  be  updated. 
To  do  this,  the  backpaths  are  iterated  again  and  for  each  node  l  encountered,  X(  )  is 
reduced  by  zf  X  CP(&). 

If  a  basis  exchange  is  required,  an  efficient  update  of  the  basis  representation  must 
preserve  the  rooted  cycle  orientation,  updating  some  entries  in  P(  )  and  D(  ).  Also, 
some  flows  X(  ),  and  some  dual  variables  U(  ),  must  be  changed.  FAC(  )  must  be 
established  for  nodes  on  newly  created  cycles.  Some  bookkeeping  in  IVAR(  )  and 
CP(  )  may  also  be  required. 

The  apparent  intricacy  of  our  task  is  deceptive.  Careful  analysis  yields  an  elegant 
solution.  However,  the  supporting  arguments  require  close  attention. 

To  simplify  the  explanation,  reorient  the  incoming  arc  ( fk,tk )  to  (i,j)  or  (/,  /)  (if 
necessary)  so  that  the  minimum  ratio  is  on  the  j  backpath.  Let  the  entering  arc  (i,j) 
have  the  outgoing  arc  ( c,d )  on  its  /'-backpath.  Also,  reorient  the  outgoing  arc  (if 
necessary)  so  that  the  first  node  encountered  on  the  j  backpath  is  c.  Figure  2  shows  a 
case  for  which  both  reorientations  are  necessary. 

In  our  example,  arc  20,  oriented  from  node  2  to  node  1 1,  leaves  the  basis  and  arc  27, 
oriented  from  node  5  to  node  13  enters.  We  call  the  backpath  segment  from  j  to  arc 
( c,d )  the  j-stem.  In  Figure  2,  i,j,  c,  and  d  are  shown.  The  /-stem  is  composed  of  nodes 
5,  2,  3,  and  11.  The  skeletal  update: 

a.  Reverses  the  orientation  of  arcs  within  the  /-stem;  and 

b.  Orients  the  entering  arc  so  that  it  precedes  this  redirected  path  with  the 
same  orientation. 

If  the  i  and  j  backpaths  merge,  the  node  where  they  merge  is  called  the  join.  The 
join  is  node  3  in  Figure  2. 

If  the  leaving  arc  lies  beyond  the  join  on  the  backpaths  a  new  cycle  is  created  in  the 
basis  exchange.  In  this  case  the  portion  of  the  i  backpath  from  node  i  to  the  node  with 
the  join  as  its  predecessor  is  called  the  i- stem.  When  node  i  is  on  the  j  backpath  the 
/-stem  is  null.  In  Figure  2,  the  /-stem  is  node  13. 

Figure  4  displays  the  pivot  logic  to  be  applied.  The  algorithm  visits  each  node 
affected  by  the  basis  update  exactly  once.  It  proceeds  up  the  /-stem  one  node  at  a  time 
visiting  the  preorder  successors  of  each  stem  node  via  IT( ).  If  a  join  is  encountered  it 
switches  to  proceed  up  the  /-stem  one  node  at  a  time  and  then  returns  to  the  /-stem.  At 
each  /-stem  node,  the  successors  of  the  next  lower  stem  node  have  already  been  visited. 
The  un visited  successors  of  the  current  stem  node  can  be  divided  into  two  groups:  the 
left  successors  are  the  nodes  visited  in  preorder  by  iterating  IT( )  from  the  current  stem 
node  until  the  next  lower  stem  node  is  encountered,  and  the  right  successors  of  the 
stem  node  are  the  remaining  unvisited  nodes  reached  by  further  iteration  of  IT(  ).  In 
Figure  2,  nodes  8,  14,  12,  and  7  are  right  successors  of  node  2  and  nodes  10,  13,  1,  6, 
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15,  9,  and  4  are  left  successors  of  node  3.  As  we  climb  the  /'-stem  the  traversal  IT( )  is 
modified  so  that  the  last  of  the  left  successors  (if  any)  points  to  the  first  of  the  right 
successors  and  the  last  of  the  right  successors  (if  any)  points  to  the  previous  stem  node. 
As  we  climb  the  y'-stem,  the  traversal  is  modified  so  that  the  last  of  the  left  successors 
(if  any)  points  to  the  first  of  the  right  successors  and  the  last  of  the  right  successors  (if 
any)  points  to  the  next  node  up  the  stem  (because  the  update  reverses  predecessor 
orientation  for  the  y'-stem). 

The  immediate  switch  to  the  /-stem  upon  encountering  the  join  during  the  y'-stem 
iteration  is  motivated  by  a  subtle  complication:  while  visiting  the  left  and  right 
successors  of  the  y'-stem  nodes,  the  nodes  on  the  /'-stem  and  their  successors  must  be 
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skipped  if  encountered.  Because  the  /-stem  nodes  are  successors  of  the  join,  visiting  the 
/-stem  as  soon  as  the  join  is  encountered  (if  one  exists)  on  the y'-stem  leaves  us  with  the 
preorder  successor  of  the  last  /-stem  node  visited.  This  valuable  artifact  enables  the 
subsequent  y'-stem  iteration  to  immediately  skip  all  /-stem  nodes  and  their  successors 
should  they  be  encountered.  This  is  the  key  step  preserving  an  efficient  one-pass  basis 
update.  In  Figure  2,  /-stem  node  13  and  its  successor  node  1  are  successors  of  the  join, 
node  3.  The  preorder  successor  of  the  last  /-stem  node  (called  the  preorder  link  in 
Figure  4)  is  node  6. 

A  stem  node  may  have  a  right  successor  which  is  on  the  root  cycle  (with  depth  zero). 
The  preorder  traversal  array  is  organized  so  that  all  successors  of  a  cycle  node  are 
encountered  before  the  next  cycle  node.  This  implies  that  a  cycle  being  broken  by  a 
leaving  arc  will  always  be  encountered  as  a  right  successor  of  a  stem  node.  In  Figure  2, 
the  broken  cycle  is  encountered  as  a  right  successor  of  node  2;  if  arc  (2, 1 1)  were  not 
the  leaving  arc,  then  node  1 1  would  be  a  preorder  successor  of  node  2. 

The  basic  arc  flows,  X( ),  are  changed  as  each  arc  is  visited  on  the  y-stem,  and  (if  a 
new  cycle  is  created)  on  the  /-stem.  If  no  cycle  is  created,  X(  )  is  changed  only  if  the 
minimum  ratio  is  nonzero,  and  then  only  on  the  arcs  visited  on  the  backpaths. 

During  the  update,  the  dual  variables  must  be  changed  so  that  ck  =  uf;  —  mkulk,  for 
every  basic  arc  k  oriented  from  node  tk  to  node/*.  With  incoming  arc  (/, y)  (reoriented 
as  in  Figure  2  so  that  the  outgoing  arc  is  on  the  y-stem),  this  relationship  is  retained  for 
all  nodes  except  for  those  which  the  update  changes  to  be  successors  of  /:  (i.e.  the 
nodes  on  the  y'-stem  and  their  successors). 

If  a  cycle  is  not  created,  this  update  proceeds  for  each  y'-stem  node  and  its  successors 
as  these  nodes  are  visited  in  preorder.  For  node  s,  and  associated  basic  arc  /  oriented 
from  s  to  p(s),  U(s)  =  C(/)  +  MUL(/)*U(P(.s)),  while  the  reverse  orientation  —  P(.s)  to 
s  yields  U(i)  =  (C(/)  -  U(-P(i)))/(-MUL(/)).  As  in  pure  networks,  the  preorder 
traversal  assures  that  a  value  of  the  dual  variable  of  a  predecessor  node  is  always 
determined  prior  to  its  use  by  any  immediate  successor  nodes. 

When  a  cycle  is  created,  the  dual  variable  must  be  determined  for  one  of  the  cycle 
nodes.  Then,  the  dual  variables  of  the  remaining  cycle  nodes  and  their  descendants 
can  be  found  one  at  a  time  in  preorder  traversal. 

This  key  dual  variable  is  computed  immediately  after  the  ratio  test  predicts  creation 
of  a  new  cycle  (the  leaving  arc  lies  beyond  the  join  on  the  ratio  backpaths).  Consider 
the  current  context  for  a  new  cycle: 


from  which  the  new  cycle  will  be  formed  (with  modified  predecessor  function  p( )): 
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The  corresponding  near- triangular  system  is: 


1 . «_,) 


a_„  b_ 


p+i 


«-3  ^  —2 


a-2 


=  C-„  + 


•  •  r_|) 


where  the  indexing  of  c  is  understood  to  yield  basis  arc  costs  (accomplished  by 
indexing  with  IVAR(  )). 

The  determination  of  one  term  (say,  ,)  of  the  solution  of  this  system  is  induced  by 
modified  forward  substitution  to  be  (e.g.,  Dantzig  1963):  w_,  =  u'_t  X  1  //,  where  (the 
cycle  factor)  /  is  defined  (again)  as 


-p 


/=  i  -  n 

8=  -  1 


and  u'  _  ,  is  the  corresponding  term  of  the  solution  of  the  strictly  upper  triangular 
system  (omitting  the  cyclic  coefficient  b_  ,  and  using  forward  substitution): 


bsu's- i 
as 


for  8  =  —  p  +  1,  .  .  .  ,  1. 


Note  that  computation  of  u'  i  requires  that  we  traverse  the  new  cycle  in  a  direction 
opposite  to  its  predecessor  orientation.  However,  before  the  update  creating  the  new 
cycle,  the  /-stem  exhibits  proper  orientation  for  at  least  part  of  our  work.  Thus,  we  can 
complete  the  first  portion  of  the  forward  substitution  for  u’_x  and  accumulate  the 
associated  partial  product  component  of  /  while  iterating  the y-stem  before  the  update. 

The  remainder  of  the  new  cycle  is  accessed  by  iterating  the  i-stem.  As  we  proceed  up 
the  i-stem  the  remaining  product  terms  of  /  are  accumulated,  and  the  reverse  /-stem 
path  is  stored  (e.g.,  using  U(  )  locations,  which  contain  obsolete  dual  values  to  be 
replaced  during  the  imminent  update).  Reaching  the  join,  this  stored  reverse  i-stem 
path  is  then  accessed  to  complete  computation  of  u'_ , . 

The  newly  created  cycle  in  Figure  5  is  composed  of  y-stem  nodes  5,  2,  and  3,  and 
i-stem  node  13.  The  dual  system  becomes 

—  »]3  +  «3 

—  «3  +  «2 

«2 

-0.88  u ,3 

which  yields 

k'i3  =  —  19.0142,  /  =  0.3488  (the  new  cycle  factor),  and 


=  45.60, 
=  15.47, 
-  0.74t/5  -  26.76, 

+  m5  =  20.67. 


Mi3  =  —54.5132  (the  new  dual  for  node  13). 

Thus,  when  a  new  cycle  is  to  be  formed,  the  new  cycle  nodes  must  be  visited  once 
(after  the  ratio  test  and  prior  to  the  pivotal  update). 

The  one-pass  preorder  traversal  update  can  now  proceed  as  presented  in  Figure  4. 
The  basis  representation  arrays  are  all  modified  on-the-fly  during  this  traversal.  The 
update  of  nodes  on  a  newly  created  cycle  (if  any)  includes  establishing  the  new  cycle 
factor  FAC(  ).  Changes  for  P(  ),  D(  ),  IT(  ),  IVAR(  ),  X(  ),  and  U(  )  proceed 
analogously  to  the  pure  network  case  (e.g.,  GNET  Bradley  et  al.  1977)  with  simple 
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(The  entering  column  is  italicized.) 


Figure  5.  New  Basis  before  Restoring  Near-Triangulation  (Entering  Column  27,  Arc  (5,  13)). 
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A  Preorder  Near-Triangulation. 


Figure  6.  New  Basis  Near-Triangulated. 
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Figure  7.  GENNET  Basis  Representation  Arrays  (for  Basis  in  Figure  5). 


modification  of  generators  for  X(  )  and  U(  )  to  accommodate  generalized  network 
coefficients  for  basic  arcs. 

Figure  5  shows  the  new  basis  (derived  from  the  example  in  Figure  2)  before  restoring 
near-triangulation  with  the  update.  A  new  cycle  is  to  be  created  and  the  new  cycle 
factor  and  a  dual  solution  (for  node  13  on  the  /-stem)  are  found  at  this  stage. 

Figure  6  displays  the  new  basis  restored  to  near-triangular  form.  At  this  point,  all 
basis  representation  arrays  are  updated  with  the  values  shown  in  Figure  7  (data  in 
italics  have  been  changed  by  the  update). 

GENNET  is  designed  to  exploit  intrinsic  pure  network  structure  commonly  found 
embedded  in  generalized  network  problems.  Note  that  when  this  algorithm  is  applied 
to  pure  network  problems,  it  automatically  adapts  to  a  minor  variant  of  GNET.  Of 
course,  GENNET  uses  floating  point  arithmetic  operations  which  are  intrinsically 
slower  on  most  computers  than  the  pure  additive  integer  arithmetic  of  GNET  (also, 
floating  point  arithmetic  requires  some  extra  editing  for  mantissa  truncation  errors). 

To  mitigate  this  disadvantage,  GENNET  can  test  logically  for  pure  network  arcs 
(with  unit  multipliers)  and  avoid  unnecessary  floating  point  multiplication  and  division 
operations. 

GENNET  also  employs  an  automatic  dual  basis  aggregation  refinement  (Bradley  et 
al.  1977,  p.  26  ff).  Explicit  values  for  D(  ),  U( )  and  IT( )  are  maintained  only  for  nodes 
with  successors.  An  array,  A(  ),  records  for  each  node  the  current  number  of  its 
aggregated  successors.  When  an  aggregated  node  is  encountered  in  the  priceout,  its 
dual  is  generated  from  that  of  the  immediate  predecessor  of  the  node.  When  a 
backpath  of  an  entering  arc  begins  with  an  aggregated  node,  it  is  disaggregated,  and 
when  the  leaving  arc  isolates  nodes  with  no  successors,  they  are  aggregated. 

Figure  8  shows  the  arrays  affected  by  the  dual  aggregation  scheme  for  the  basis  in 
Figures  2  and  3.  IT(  )  indicates  aggregated  nodes  with  0  entries,  and  these  nodes  have 
broken  outlines  in  the  one-tree  depiction.  The  priceout  of  arc  (5, 13)  requires  that  U(5) 
be  generated  using  the  predecessor  dual  U(2).  Node  5  subsequently  starts  the  j- 
backpath  and  is  disaggregated.  The  update  leaves  node  1 1  with  no  successors,  and  thus 
aggregated. 

5.  Computational  Experience 

Significant  design  alternatives  for  GENNET  have  each  been  evaluated  by  extensive 
experimentation  at  large  scale.  Illustrative  computational  experience  is  abstracted  in 
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0 
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0 

3 
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5 
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0 

13 
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4 
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0 
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0 

9 
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11 
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0 

Aggregated  One-Tree 

Figure  8.  Aggregated  Basis  Representation  (for  Figures  2  and  3). 

this  section  for  some  of  the  prototypic  systems  tested.  Departing  somewhat  from  the 
style  of  the  paper  documenting  such  work  for  GNET  (Bradley  et  al.  1977),  relative 
performance  is  reported  even  for  some  competitive  design  features  subsequently 
rejected  for  adoption  (some  readers  of  Bradley  et  al.  1977  have  concluded,  quite 
incorrectly,  that  only  those  features  reported  for  GNET  were  tested).  This  should  help 
other  researchers  avoid  our  mistakes,  and  may  even  change  some  widely  held  miscon¬ 
ceptions  and  correct  a  few  translation  errors  in  textbooks. 

Among  the  key  issues  to  be  resolved  are: 

(a)  Static  Storage  of  Arcs.  The  arc  lists  can  be  stored  in  arbitrary  sequence,  or,  to 
save  space,  arcs  can  be  stored  contiguously  by  tail  node,  or  by  head  node,  thereby 
replacing  an  arc-length  index  array  by  a  (presumably  much  shorter)  hierarchical 
node-length  entry  point  array. 

(b )  Preorder  Manipulation  of  Basis.  The  triple  label  (augmented  predecessor  index) 
method  (Glover  et  al.  1972,  1973)  will  be  presented  and  compared  with  the  preorder 
traversal  method. 

(c)  Basis  Aggregation.  An  aggregated  basis  representation  will  compete  with  an 
explicit  representation. 

(d)  Pricing  Schemes.  Candidate  list  schemes  and  explicit  arc  pricing  mechanisms 
widely  used  in  general  linear  programming  systems  will  vie  with  dynamic  candidate 
queue  disciplines. 
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200 
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60 
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8 

60 
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80 
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12 
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80 
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50 
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50 

50 

3,400 

0 

NG30 

1,000 

50 

50 

4,400 

0 
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50 
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0 

NG32 

1,500 

75 

75 

4,342 

0 

NG33 

1,500 

75 

75 

4,385 

0 

NG34 

1,500 

75 

75 

5,107 

0 

NG35 

1,500 

75 

75 

5,730 

0 

NETGENG 

(Generalized) 

GT01 
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100 

100 

1,500 

0 

GT02 

200 

100 

100 

2,000 

100 

GT07 

300 

135 

115 

4,000 

0 

GT12 

400 

20 

100 

5,000 

0 

GT15 

1,000 

5 

995 

4,000 

100 

GT16 

1,000 

20 

100 

6,000 

100 

GT18 

1,000 

30 

400 

7,000 

0 

Figure  9.  Some  Benchmark  Problems  (Problem  NG27  is  extensively  studied  in  Bradley  et  al.  1977). 


(e)  Pure  Network  Specialization.  Generalized  network  algorithms  would  ideally  adapt 
to  pure  networks  with  efficiency  comparable  to  pure  network  codes. 

(f)  Starting,  Tuning  and  Tailoring.  Which  algorithm  parameters  and  settings  lead  to 
high  efficiency  for  interesting  classes  of  problems?  Are  heuristics  for  advanced  starting 
solutions  worthwhile? 

(g)  Generalizations.  Advanced  features  and  generalizations  will  be  suggested. 

Computational  tests  have  been  made  with  many  problems,  including  a  benchmark 

suite  of  pure  network  problems  generated  with  NETGEN  (Klingman  et  al.  1974),  and 
generalized  network  problems  generated  by  NETGENG  (Glover  et  al.  1977,  1978). 
Figure  9  gives  some  problem  characteristics. 

Static  arc  storage  has  been  implemented  in  three  ways: 

•  Contiguous  by  head  node  with  hierarchical  node-length  entry  point  array, 

•  Contiguous  by  tail  node,  with  hierarchical  node-length  entry  point  array,  and 

•  Explicit  arcs  in  arbitrary  order. 

In  competition  with  our  basis  manipulation  using  preorder  traversal,  a  triple  label 
representation  originally  suggested  by  E.  Johnson  (1966)  has  been  implemented  for 
pure  networks  and  called  the  augmented  predecessor  indexing  method  by  Glover, 
Karney  and  Klingman  (1972).  Glover,  Klingman  and  Stutz  (1973)  report  that  the 
method  has  been  extended  to  generalized  networks,  but  reveal  no  details. 

We  have  implemented  our  own  efficient  version  of  the  triple  label  scheme. 

The  triple  label  representation  uses  predecessor,  successor,  and  brother  pointers  for 
each  node.  Figure  10  shows  these  arrays  for  the  basis  in  Figure  2. 

To  briefly  illustrate  the  triple  label  scheme  for  generalized  networks,  let  our  situation 
be  the  same  as  for  the  preorder  example  (the  y'-backpath  of  the  incoming  arc  is 
/  arranged  to  include  the  outgoing  arc  and  to  encounter  node  c  first  on  that  backpath). 
Let  vd+ ,  =  i,  and  the  backpath  from  y  to  c,  vd,vd_v  .  .  .  ,  c.  The  skeletal  update 
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Figure  10.  Triple  Label  Basis  Representation  (for  Figure  2). 


scheme  is: 

1.  Set  S  =d. 

2-  If  S(o4_,)  =  vs,  go  to  Step  3. 

Otherwise,  if  possible,  find  a  node  v*  such  that  P(o*)  =  vs_  {  and  B(o*)  =  v$,  and  set 
B(t>*)  <—  B(t>8),  Go  to  Step  4. 

3.  Set  S(«g_i)<-B(u5). 

4.  Set  P(«8)<-t>8+1,  B(o8)<-S(c8+1),  S(%+  [)<-  o8. 

5.  Set  8<-  8  —  1.  Then,  if  vs+ ,  ^  c,  Go  to  Step  2.  Otherwise,  Stop. 

(Step  2  exhibits  the  key  extension  for  generalized  networks  of  the  pure  network  triple 
label  scheme.) 

These  five  steps  reverse  the  orientation  of  all  arcs  on  the  y'-stem  and  orient  (vd+ , ,  vd) 
so  that  it  begins  this  redirected  path.  All  other  triple  label  operations  are  obvious 
alterations  of  the  preorder  traversal  procedure. 

A  static  candidate  list  pricing  strategy  (e.g.,  Glover  et  al.  1977,  1978,  Mulvey  1978) 
an  explicit  arc  pricing  method  reminiscent  of  general  linear  programming  systems,  and 
a  dynamic  candidate  queue  (Bradley  et  al.  1977)  have  been  tested. 

The  ( Lx,L2 )  candidate  list  procedure  is  a  simple  strategy.  Nodes  with  leaving  arcs 
(or  entering  arcs  for  contiguous  head  node  arc  lists)  are  sequentially  priced,  placing  the 
most  negative  candidate  (if  any)  from  each  node  on  the  candidate  list  until  L2 
candidates  have  been  located  (arbitrarily  organized  arc  lists  require  explicit  arc 
pricing).  Entering  arcs  are  chosen  from  the  candidate  list  by  most  negative  reduced 
cost  until  L,  iterations  have  been  carried  out,  or  until  all  list  entries  have  nonnegative 
reduced  costs.  Arcs  with  nonnegative  reduced  costs  are  dropped  from  the  list  and 
replaced  by  continuing  to  scan  the  nodes  until  L2  candidates  have  been  found.  If  an 
exhaustive  pass  through  the  nodes  results  in  less  than  L2  candidates,  then  an  optional 
closing  gambit  sets  L2  equal  to  the  actual  candidates  found,  and  reduces  L,  by  half 
unless  L,  equals  one.  Glover  et  al.  (1977,  1978)  report  (L,,L2)  of  (5, 10)  to  be  best  in 
their  work. 

The  candidate  queue  is  a  dynamic  list  of  interesting  arcs  and  nodes,  scanned  in  a 
cyclic  manner.  The  entering  arc  is  selected  from  the  queue  by  pricing  NNE  entries;  if 
an  interesting  node  is  encountered  it  is  replaced  by  its  best-priced  entering  arc  (or 
leaving  arc  for  contiguous  tail  node  arc  lists).  Arcs  pricing  favorably  are  retained  in  the 
queue.  When  the  end  of  queue  is  encountered,  the  queue  is  refreshed  by  pricing  IPG 
nodes  in  a  cyclic  general  arc  scan.  During  an  opening  gambit  of  NNS  pivots,  the  nodes 
incident  to  the  entering  basic  arc  are  added  to  the  queue.  There  is  no  closing  gambit, 


oo 


Generalized  Network  Codes 

Pure  Network  Codes 

(start) 

GNET 

GNET 

Problem 

(5, 10) 

NNE  =  32 

NNE  =  32 

NNE  =  32 

NNE  =  8 

NNE  =  16 

NNE  =  8 

NETGEN 

TLA 

TQA 

HQA 

HQP 

HQPX 

(Big-M) 

HQP 

HQPX 

SUPERK 

NG15 

2.19 

2.17 

3.10 

2.50 

2.07 

1.5 

1.37 

1.17 

2.15 

NG18 

0.83 

0.60 

0.63 

0.63 

0.44 

2.5 

0.42 

0.39 

1.30 

NG19 

1.71 

1.10 

0.90 

0.79 

0.70 

2.0 

0.67 

0.56 

2.36 

NG22 

0.97 

0.61 

0.57 

0.64 

0.46 

2.0 

0.43 

0.38 

1.54 

NG23 

1.60 

0.76 

0.84 

0.63 

0.71 

1.5 

0.48 

0.43 

2.11 

NG26 

1.03 

0.48 

0.54 

0.52 

0.45 

2.0 

0.36 

0.30 

0.93 

NG27 

1.11 

0.74 

0.78 

0.71 

0.59 

1.5 

0.48 

0.47 

1.17 

NG28 

4.42 

2.20 

2.71 

2.13 

1.67 

2.5 

1.41 

1.36 

4.53 

NG29 

5.96 

2.32 

3.12 

1.96 

1.95 

2.0 

1.51 

1.40 

4.42 

NG30 

4.52 

2.58 

3.00 

2.27 

2.30 

1.5 

1.67 

1.45 

5.58 

NG31 

5.46 

2.58 

3.40 

2.27 

2.53 

1.5 

1.69 

1.58 

5.47 

NG32 

9.79 

3.32 

4.42 

3.32 

3.12 

2.0 

2.56 

2.28 

7.92 

NG33 

9.89 

4.18 

4.46 

3.15 

3.49 

2.5 

2.61 

2.58 

6.83 

NG34 

9.22 

4.24 

4.83 

3.47 

3.53 

1.5 

2.94 

2.44 

8.94 

NG35 

10.18 

3.83 

5.67 

3.59 

3.90 

1.5 

2.77 

3.00 

9.32 

Seconds: 

68.88 

31.71 

38.97 

28.58 

27.91 

21.37 

19.79 

64.57 

Rank: 

5 

3 

4 

2 

1 

Legend:  Head/Tail  arc  organization 

List/Queue  pricing  mechanism 
API  triple  label /Preorder  traversal 
X  aggregated  successors 


Figure  11.  Pure  Network.  Test  Problems. 
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6.41 
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GT07 

2.30 

2.91 

3.03 

0.78 

1.63 

0.94 

0.82 

0.56  ! 

2 

GT12 

3.62 

4.59 

3.34 

4.01 

3.66 

3.49 

3.02 

2.77  ' 

1.5 

GT15 

13.41 

24.14 

28.99 

6.46 

9.92 

8.43 

7.94 
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GT16 

9.98 

17.84 

7.47 

12.08 

12.16 

9.84 

5.77 
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GT18 

12.71 

19.19 

13.63 

17.89 

13.89 

13.05 

10.46 

9.72  j 
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Seconds: 

47.46 

74.44 

62.34 

50.36 

46.39 

41.67 

34.18 

28.15 

Rank: 

5 

8 

7 

6 

4 

3 

2 

1 

Legend:  Head/Tail/Explicit  arc  organization 
List/Queue  pricing  mechanism 
API  triple  label/ Preorder  traversal 
X  aggregated  successors 

Figure  12.  Generalized  Transshipment  (GT)  Test  Problems. 
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since  the  queue  automatically  shrinks  and  finally  collapses  at  optimality.  Bradley, 
Brown,  and  Graves  (1977)  suggest  NNE  =  32,  NNS  =  3m/ 4  and  IPG  =  m/10  +  1  for 
pure  network  problems  with  m  nodes. 

A  rule  to  break  ties  in  the  ratio  test  which  guarantees  finite  convergence  for  pure 
transshipment  problems  has  been  developed  by  Cunningham  (1976).  Bradley,  Brown 
and  Graves  (1977)  show  that  the  conditions  necessary  for  finite  convergence  are 
naturally  satisfied  by  GNET  on  over  90%  of  its  degenerate  pivots.  Elam,  Glover  and 
Klingman  (1979)  have  observed  that  the  results  of  Cunningham  can  only  be  extended 
to  the  generalized  network  case  when  the  multipliers  are  positive.  We  have  not  used 
Cunningham’s  modification. 

Bland  (1977)  presents  a  class  of  restrictions  of  pricing  and  ratio  tests  for  general 
linear  programs  which  relies  exclusively  on  primal  simplex  representation  and  guaran¬ 
tees  finite  convergence.  These  rules  are  easily  modified  to  produce  an  efficient  finite 
simplex  algorithm.  The  modification  interferes  with  effective  pricing  strategies  only 
during  degenerate  pivot  sequences,  and  the  restrictions  increase  in  severity  only  with 
the  number  of  pivots  in  that  sequence.  However,  during  a  degenerate  pivot  sequence 
restriction  records  must  be  accumulated  (e.g.,  a  list  with  each  incoming  variable  in  one 
of  our  schemes).  This  record  is  naturally  accommodated  by  the  dynamic  candidate 
queue,  but  not  by  a  static  candidate  list  or  explicit  pricing.  No  purpose  is  served  by 
reporting  such  unbalanced  competition. 

A  starting  strategy  has  also  been  tested  in  conjunction  with  pricing  alternatives.  A 
straightforward  starting  method  examines  each  node  with  supply  (or  demand  for 
contiguous  arc  storage  by  head  node)  and  assigns  as  much  flow  as  possible  to  its  least 
cost  leaving  (entering)  arc.  The  procedure  stops  when  an  exhaustive  pass  of  the  nodes 
makes  no  additional  flow  assignments.  The  starting  solution  achieved  is  not  necessarily 
feasible  (e.g.,  “exhaustive  pass  sequential  source  minimum  start”  Glover  et  al.  1978). 

Artificial  arcs  are  driven  from  the  basis  in  all  experiments  using  a  Big-M  method 
(e.g.,  Bradley  et  al.  1977).  This  choice  is  principally  motivated  by  the  comparability  of 
competitive  tests  between  pure  and  generalized  network  codes  on  pure  and  generalized 
network  problems.  (A  two-phase  method  is  employed  in  production  use.) 

Choosing  the  best  Big-M  value  is  a  bit  tricky.  The  smallest  Big-M  value  which  yields 
a  feasible  optimal  solution  (if  one  exists)  is  best  in  our  experience.  Small  Big-M  values 
may  fail  to  produce  feasibility,  and  large  values  inflict  numerical  difficulties.  In 
practice,  a  default  value  is  used  and  an  automatic  restart  recovery  is  applied  if  an 
infeasible  solution  persists.  If  a  restart  with  a  higher  Big-M  value  fails  to  reduce  the 
total  infeasibility,  a  terminally  infeasible  solution  is  declared.  Figures  11  and  12 
indicate  the  multiple  of  maximum  absolute  arc  cost  used  for  Big-M  in  each  problem. 

Computational  tests  have  been  performed  on  various  computer  systems.  The  times 
reported  here  are  accurate  to  the  precision  displayed  for  IBM  370/168-3  using  the 
FORTRAN  H  compiler  with  OPT(2).  Solution  times  exclude  input/output  overhead. 
GENNET  uses  double  precision  (IBM  REAL*  8)  arithmetic  for  floating  point  opera¬ 
tions  and  storage. 

Solution  times  are  given  in  Figure  1 1  for  the  pure  network  test  problems.  Perfor¬ 
mance  is  given  for  three  pure  network  codes  (two  versions  of  GNET  and  SUPERK)  as 
well  as  for  several  representative  generalized  network  prototypic  systems.  We  can  thus 
compare  the  best  generalized  network  scheme  with  the  best  pure  network  code 
exhibiting  equivalent  features  (GNET,  or  its  aggregated  successor  variant  (Bradley  et 
al.  1977). 

The  times  for  SUPERK  also  provide  the  only  available  objective  means  for 
comparison  of  our  implementations  to  that  of  Glover,  Hultz,  Klingman  and  Stutz 
(1973);  they  report  that  their  generalized  network  code,  NETG,  is  about  as  fast  as 
SUPERK  (a  fast  out-of-kilter  code  for  pure  networks)  when  both  are  used  to  solve 
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Problem 

Nodes 

Arcs 

Seconds 

Pivots 

AIRLP 

170 

3,040 

2.62 

420 

COAL 

170 

3,923 

1.80 

471 

STEEL 

422 

1,279 

0.39 

499 

FOAM 

951 

4,953 

3.74 

1,258 

ODSAS(GN) 

1,431 

4,615 

3.22 

1,427 

ALUMINUM(GN) 

2,178 

7,216 

3.57 

2,794 

REFINE(GN) 

3,110 

6,617 

4.72 

3,322 

FOOD(GN) 

3,716 

13,907 

12.11 

7,004 

(Big-M  =  10  x  largest  cost  coefficient) 


Figure  13.  (GN)  Test  Problems  (GN  rows  extracted  from  real-world  LP/MIP  models  Brown,  McBride 
and  Wood  to  appear  with  null  columns  deleted  and  slack  arcs  added). 

pure  networks.  Using  their  version  of  SUPERK  on  our  computer  we  have  shown  that 
our  implementation  of  NETG  (called  TLA  in  our  nomenclature)  is  at  least  as  efficient 
as  their  claims  for  NETG. 

However,  in  these  tests  TLA  is  substantially  outperformed  by  the  alternate  systems. 
For  the  pure  network  problems,  best  performance  is  achieved  by  (Figure  11): 

Contiguous  arcs  by  head  node, 

Candidate  queue  pricing, 

Preorder  traversal,  and 

Aggregated  successors. 

TLA  does  not  incorporate  any  of  these  features,  and  is  generally  less  than  half  as  fast 
as  competitors. 

Although  GENNET  (HQPX)  should  in  theory  rival  GNET  (HQPX)  with  pure 
network  problems,  the  overhead  of  testing  in  GENNET  for  more  general  basis 
structure  and  the  additional  computational  burden  of  floating  point  arithmetic  exact  a 
performance  penalty  of  about  15  percent. 

Figure  12  shows  solution  times  for  the  generalized  transshipment  network  problems. 
Note  that  the  starting  strategy  helps  candidate  list  performance  and  hinders  the 
candidate  queue.  Arranging  arcs  contiguously  by  head  node  dominates  both  tail  node 
and  explicit  arc  list  designs.  The  candidate  queue  provides  good  performance  if 
accompanied  with  contiguous  arcs  by  head  node.  Preorder  traversal  continues  to 
provide  better  performance  than  triple  label  representation  in  all  design  contexts. 
Aggregated  successors  offer  a  pronounced  advantage. 

GENNET  (HQPX)  provides  best  overall  performance.  It  offers  a  decided  advantage 
on  problems  with  many  more  sinks  than  sources,  a  situation  common  in  real  life. 

Figure  13  displays  performance  of  GENNET  (GN,  HPQX)  applied  to  a  set  of  (GN) 
problems  extracted  from  a  collection  of  real-world  LP/MIP  models  (Brown,  McBride 
and  Wood  to  appear).  Despite  the  slight  additional  floating  point  arithmetic  required 
to  solve  (GN),  GENNET  solves  these  problems  much  faster  than  would  be  predicted 
by  experience  with  randomly  generated  (GT)  problems.  Tuning  of  the  pricing  mecha¬ 
nism  greatly  enhances  this  difference. 

Close  scrutiny  of  solution  trajectories  lends  some  insight  into  GENNET’s  good 
performance.  GENNET  has  a  one-pass  iteration  unless  a  cycle  is  formed;  a  cycle  is 
formed  on  only  5-to-24  percent  of  all  iterations  for  these  problem  sets — 5-to-10  percent 
for  most  problems.  Also,  the  explicit  (nonaggregated)  subset  of  the  nodes  is  remark¬ 
ably  small,  seldom  numbering  much  more  than  the  number  of  source  nodes.  Finally, 
the  length  of  backpaths  is  quite  short,  averaging  about  the  number  of  echelons  (path 
length  from  sources  to  sinks)  in  the  model,  or  just  more  than  2  in  these  problems. 

6.  Conclusion 

The  generalized  network  system  GENNET  is  small,  fast  and  easy  to  modify.  Adaptations  have  already 
included  using  GENNET  in  a  system  to  solve  generalized  networks  with  complicating  side  constraints 
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and/or  complicating  variables  (McBride  to  appear).  GENNET  has  also  been  incorporated  in  a  powerful 
microcomputer-based  network  optimization  system  by  Brown,  Duff  and  Finley  (1982)  and  Finley  (1982) 
using  an  APPLE-II  host  and  PASCAL  implementation  language.  Modifications  for  mixed  integer  general¬ 
ized  networks  have  also  been  tested  (though  not  with  care  sufficient  to  warrant  publication  at  this  time). 
GENNET  has  proven  to  be  a  worthy  successor  of  GNET  (Bradley  et  al.  1977). 

Preorder  traversal  is  appealing  for  its  mathematical  and  implementation  elegance,  and  has  proven  to  be 
efficient  and  flexible  for  generalized  networks  (as  it  was  for  pure  networks).  (Adolphson  and  Heum  1981 
have  also  suspected  this  and  have  independently  pursued  this  avenue.) 

Experience  shows  that  the  GENNET  design  performs  much  more  efficiently  on  real  models  than  on 
randomly  generated  test  problems  of  nominally  equivalent  size;  this  design  is  also  technically  and  philosoph¬ 
ically  compatible  with  the  various  systems  we  have  devised  for  solving  other  more  general  classes  of 
optimization  models. 

The  FORTRAN  programs  GENNET — (GT)  and  (GN)  versions — (Copyright  1984)  are  licensed  to 
researchers  for  a  nominal  charge  on  an  exclusive  use  basis.  For  further  information  write  the  authors  via 
P.O.  Box  1832,  Alexandria,  Virginia  22314,  USA. 
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